Determining by Raman spectroscopy the average thickness and N-layer-specific surface coverages of MoS2 thin films with domains much smaller than the laser spot size

Raman spectroscopy is a widely used technique to characterize nanomaterials because of its convenience, non-destructiveness, and sensitivity to materials change. The primary purpose of this work is to determine via Raman spectroscopy the average thickness of MoS2 thin films synthesized by direct liquid injection pulsed-pressure chemical vapor deposition (DLI-PP-CVD). Such samples are constituted of nanoflakes (with a lateral size of typically 50 nm, i.e., well below the laser spot size), with possibly a distribution of thicknesses and twist angles between stacked layers. As an essential preliminary, we first reassess the applicability of different Raman criteria to determine the thicknesses (or layer number, N) of MoS2 flakes from measurements performed on reference samples, namely well-characterized mechanically exfoliated or standard chemical vapor deposition MoS2 large flakes deposited on 90 ± 6 nm SiO2 on Si substrates. Then, we discuss the applicability of the same criteria for significantly different DLI-PP-CVD MoS2 samples with average thicknesses ranging from sub-monolayer up to three layers. Finally, an original procedure based on the measurement of the intensity of the layer breathing modes is proposed to evaluate the surface coverage for each N (i.e., the ratio between the surface covered by exactly N layers and the total surface) in DLI-PP-CVD MoS2 samples.


Other intensity references
In the main text, the MoS2 A1g and E 1 2g modes integrated intensities are normalized relatively to the integrated intensity of the 521 cm -1 mode of a bare Si(111) wafer with only native oxide.In the literature or in future works other intensity references can be considered and Raman set-ups with different characteristics can be used.To enable results comparison, we have measured i) bulk MoS2 crystal (HQ graphene) for which we found, A(A1g)/A(Si111) = 0.116 and A(E 1 2g)/A(Si111) = 0.067, ii) the 521 cm -1 mode of Si(100) with 90 ± 6 nm SiO2 (resp.native oxide) which integrated intensity equals 1.32 ± 0.05 (resp.0.68) relative to Si(111) and iii) the polarization ratio of our collection and detection system which is H:V=1:0.76,with H (resp. V) standing for horizontal (resp.vertical) which is parallel (resp.perpendicular) to laser polarization.S4

2D growth toy model
Contrary to standard CVD growth of MoS2, the DLI samples were produced rather quickly with a large amount of precursors, and the samples have small domain sizes (around 50 nm).The model is thus based on a simplifying assumptions 1 : 1) Once a flake has nucleated, the advance rate of the growth front 2   = 1 nm/s is constant (only dictated by the incorporation to the flake edge).When another flake is met the advance stops.
The "high-performance grid-based spatial simulation framework" DynamicGrids.jl[3,4] was chosen for its speed and simplicity.It is part of the Dispersal.jlframework [3].At the time of writing only square grids were readily available, but the goal was only to show why multilayers seemed to "wait for the previous layer to grow" before appearing.The model shows that it is possible to qualitatively reproduce these observations with constant probabilities.The resulting shapes are squared instead of hexagonal or triangular, but that changes only details such as the relation between characteristic size and perimeter or 1 Ref. [2] also used a constant advance rate, surface diffusion seemed fast enough to be ignored even for micrometer-sized flakes in the conditions they considered. 2Ref. [2] Kinetic Monte Carlo model also showed a "line-by-line" growth when the flake is large enough for the atomistic irregularity of the edge to be indiscernible.In our case, it is the "collisions" with other nano-domains that lead to irregular shapes and smooth out the curves.S7 area.Those details would be important if the goal was to extract quantitative values of the parameters, which is out of scope here (much more elaborate models exist, see e.g.[2]).
The sample area is divided in   ×   cells on a square grid.Each cell holds the number of layers covering it.So at each time step (time ), the sample is described by a   ×   matrix of integers.
The matrix is updated for the next time step (time  + ) as follow, for each cell: 1) If at time  any of the cell first neighbors is higher than the cell (meaning there is a growth edge nearby), then increment the cell.So the nearby edge advances to the cell.Note: it does not matter if there are several edges nearby; the cell is incremented only once.
2) Otherwise, randomly increment the cell by 1, with a fixed probability   =  2 .This crudely models the germination process.
The time increment  = 1 s, and the cell size  = 1 nm define the linear advance rate   = /.For a 2 min process (120 s) on a 2000x2000 grid the calculation takes less than a second on an AMD Ryzen 9 3900X 12-Core processor.A recording of the growth is available as "growth.mp4".
Of course, there are some run to run fluctuations in the coverage evolution, which would be akin to fluctuations from one sample position to another.But fluctuations happen to be negligible in this case, and the grid area (2 μm square) is still smaller than the probed area (the spectra are averaged over several spot sizes to increase signal to noise ratio).This would have to be reconsidered if the spatial resolution were not limited by diffraction (with a tip-enhanced Raman spectroscopy probe for instance), or if the domain sizes were larger (fewer domains probed at once).

S8
This set of parameters seemed reasonable and gave an acceptable agreement with the experiments.
The curves in the following figure are remarkably robust to any parameters change (neither doubling or halving the cell size, and thus the advance rate, nor multiplying or dividing by 5 the growth rate).
In the following, the capital  stands for number of layers, while the lower case   stands for the "number of cells that are covered by exactly  layers".We'll call   the total number of cells used in the calculation.With these definitions, the coverage by exactly  layers is and the average number of layers is 〈〉 =

Figure S1 : 2 .
Figure S1: (a) Optical microscopy picture of a mechanically exfoliated MoS2 sample on a

Figure S4 :
Figure S4: Representative AFM images of DLI-PP-CVD samples with  ̅ < 1.25 ((a) and (d)), 1.25 <  ̅ < 1.3 ((b) and (e)) and 1.3 <  ̅ < 2 ((c) and (f)).(a), (b) and (c) topological are scarce, the  ≥ 2 terms are negligible and 〈〉~ 1 .This explains why close to 0 the  1 (〈〉) slope is 1, and  0 (〈〉)~1 − 〈〉.Two effects tend to reduce the slope of  1 (〈〉) : (i) the growth of the second layer reduces  1 , since a cell with two layers is no longer counted as a monolayer.(ii) at some point a monolayer flake encounters other monolayers, which limits their expansion.This is why monolayers eventually get covered entirely ( 1 → 0), despite the fact that all layers have the same advance rate in our simplified model.This model also provides some support to the Raman scattering analysis made in the article core as shown in FigureS5.For 〈〉 ∈ [1.25, 1.3] the  = 1 and  = 2 contributions indeed dominate.For 〈〉 ∈ [1.3, 2] the  = 3 contribution rises.For 〈〉 ∈ [2, 2.75] the  = 1 or the  = 4 contributions are to be taken into account, in addition to the  = 2 and  = 3.This most challenging region is thus treated last in the proposed analysis.Finally, for 〈〉 ∈ [2.75, 2.85] both the modeled  = 1 and  = 5 contributions are negligible.S9 So the assumptions underlying the proposed procedure are valid in the model.

Figure S5 :
Figure S5: Coverage   (ratio of the surface covered by exactly N layers to the total surface) as a function of the average number of layers.The vertical dashed lines mark the